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1. LITHOSPHERIC STRUCTURE BENEATH ARABIA 

T. A. Mokhtar, C. J. Ammon, R. B. Herrmann, and H. A. A. Ghalib 

1.1 Abstract 

The three-dimensional seismic velocity distribution beneath the Arabian 
Plate is investigated using Love and Rayleigh waves. We obtained a balanced 
path coverage using seismograms generated by earthquakes located along the 
plate boundaries in the Red Sea, Gulf of Aqaba, Gulf of Aden, western Iran, 
Turkey, and the Dead Sea fault system. We measured Love- and Rayleigh- 
wave group velocity dispersion using multiple-filter analysis and then per¬ 
formed a tomographic inversion using these observations to estimate lateral 
group-velocity variations in the period range of 5-60 seconds. The Love- and 
Rayleigh-wave results are consistent and show that the average group veloc¬ 
ity increases from 2.38 and 2.44 km/ a at 5-7 seconds to 3.74 and 3.98 at 56-60 
seconds, for Rayleigh and Love waves respectively. The tomographic results 
also delineate first-order regional structure heterogeneity as well as a sharp 
transition between the two major tectonic provinces in the region, the Ara¬ 
bian Shield (faster than average) and the Arabian Platform (lower than aver¬ 
age). 

1.2 Introduction 

The deployment of broadband seismic stations within the Arabian Shield 
(Vernon and Berger, 1997) provided an excellent opportunity to study the 
seismic structure of the Arabian plate using high-quality seismic signals that 
were not previously available for this part of the world. A number of recent 
studies have made use of the recorded broadband data (e.g. Sandvol et al., 
1998, Mellors et al., 1997, Baker et al., 1997, McNamara et al., 1997, Rodgers 
et al., 1997, Mokhtar et al., 1997). In this paper we present maps of the aver¬ 
age group velocity variation across the Arabian Plate. One goal of our effort 
was to increase the lateral resolution by incorporating as many as possible of 
the available wave propagation paths across the Arabian plate, and espe¬ 
cially, incorporate observations in which the records were made within the 
Arabian plate. In the period range between 5 and 60 seconds we had a maxi¬ 
mum of 916 Rayleigh and 653 Love waves. Since we relied on short paths and 
small events, the fewest paths are those at the longest periods where we col¬ 
lected 376 Rayleigh wave and 178 Love- wave observations. 

1.3 Tectonic Setting 

The Arabian Plate consists of two major tectonic provinces, the Arabian 
shield and the Arabian Platform (Figure 1). The Arabian shield covers about 
one third of the Arabian Peninsula and consists of Precambrian gneiss and 
metamorphosed sedimentary and volcanic rocks that have been intruded by 
granites (Powers et al., 1966, Brown, 1972). The shield consists of five micro¬ 
plates (Afif, AR Rayn, Asir, Midyan, and Hijaz micro-plates (Stoeser and 
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Camp (1985)) which are separated by four ophiolite-bearing suture zones. 
These micro-plates are considered the remnants of Precambrian island arcs 
(Schmidt et al., 1979) that accreted to form an Arabian neocraton around 
630 Ma and which in turn were subjected to subsequent intra-cratonic defor¬ 
mation and magmatism producing the present day shield (Stoeser and Camp 
(1985)). Widespread Tertiary and Quaternary volcanic rocks related to initial 
stages of the Red Sea fQrmation are predominant along western Arabia in the 
shield (Brown, 1972, Coleman, 1977). 
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Figure 1. Tectonic provinces of the Arabian plate and the different plate boundaries. 

The Arabian platform is a large sedimentary basin that comprises about 
two thirds of the plate and consists of Palaeozoic and Mesozoic sedimentary 
layers that uncomformably overlap the basement rocks and gently dip to the 
east (Powers et al., 1966). Platform sediments increase in thickness to the 
east and reach a thickness of 10 km or more beneath the Mesopotamian fore¬ 
deep (Brown, 1972). The western platform sediments are relatively 
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undeformed but deformation increases to the east, towards the foredeep and 
the Zagros and Taurus mountains. Five tectonic boundaries surround the 
Arabian plate (Maamoun, 1976): The continental collision boundary between 
the Arabian plate and the Persian and Turkish plates along the Zagros and 
Taurus mo untains in the north and northeast; the subduction boundary in 
the Gulf of Oman region in southern Iran (Makran subduction zone); the 
transform fault boundary along the Owen Fracture zone in the southeast; 
the transform fault boundary along the Dead Sea fault system in the north¬ 
west; and the spreading axis along the Red Sea and Gulf of Aden in the west 
and southwest. All of these boundaries are tectonically active and produce an 
appreciable number of earthquakes, especially along the boundary associated 
with the Zagros mountain belt in the northeast. 

1.4 Previous Geophysical Investigations 

The shear-wave velocity structure of the two major tectonic provinces 
was modeled by Mokhtar & Al-Saeed (1994) using a surface-wave dispersion 
inversion. The model of the Arabian platform was found to be similar to that 
of East Africa and consists basically of an upper crust, 20 km thick, with 
shear wave velocity of 3.4 km/s overlying a 20 km thick lower crust with 
shear wave velocity of 4.0 km/s. Mokhtar (1995) used waveform modeling to 
verify the Arabian platform shear-wave velocity models of Mokhtar and Al- 
Saeed (1994). The Arabian shield velocity model consists of an upper and 
lower crust of comparable thickness to those of the platform with P-wave 
velocities of 6.3 km/s and 7.0 km/s for the upper and lower crust respectively 
(Mooney et al., 1985) and S-wave velocities of 3.6 km/s and 3.88 km/s for the 
upper and lower crust respectively (Mokhtar and Al-Saeed (1994)). The depth 
to the mantle is about 45 km beneath the platform and decreases to the 
southwest reaching about 40-38 km in the southwestern part of the shield 
(Mokhtar and Al-Saeed (1994), Mooney et al., (1985), Badri (1991)). Sandvol 
et al., (1998) estimated the lithospheric mantle and crustal velocity struc¬ 
tures beneath the Arabian shield through the modeling of teleseismic P 
waves recorded by the temporary broadband array used in this study. Appli¬ 
cation of the receiver function techniques showed that the crustal thickness 
of the shield area varies from 35 to 40 km in the west adjacent to the Red 
Sea, to 45 km in central Arabia. These results are consistent with previous 
results from surface-wave inversion (Mokhtar et al., (1994), and the deep seis¬ 
mic refraction results (Mooney et al., (1985) and Badri (1991)). Ghalib et al., 
(1998) and Ghalib (1992) utilized Rayleigh-wave fundamental-mode group- 
velocity observations from five analog seismic stations to investigate the 
three-dimensional seismic structure of the Arabian plate. They reported the 
presence of two discontinuities at 15-22 and 35-55 km depth and the crustal 
velocity was found to be higher under the Arabian shield than the rest of the 
plate. Ritzwoller and Levshin (1998) and Ritzwoller et al., (1998) produced 
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tomographic maps from surface-wave group velocities across all of Eurasia. 
These maps were at a length scale intermediate between regional and global 
surface-wave studies. 

Mokhtar and Al-Saeed (1994) presented two sets of dispersion data (Path 
II and Path III) for Love and Rayleigh-waves for the Arabian platform and 
one set for Rayleigh-waves for the Arabian shield (Path I). Path II extends 
from the eastern part of the Gulf of Aden to RYD station, while Path III con¬ 
nects the events in southern Iran to RYD station. There were no reliable 
Love-wave dispersion data for the Arabian shield reported by Mokhtar and 
Al-Saeed (1994). The Arabian shield has a higher group velocity than the 
average while the Arabian plate group velocity values are lower than the 
average. Similar behavior is observed for Love-waves. 

1.5 Observations 

We measured surface-wave group velocities generated by earthquakes 
located along the boundaries of the Arabian plate in the Red Sea, Gulf of 
Aqaba, Gulf of Aden, western Iran, Turkey, and the Dead Sea fault system. 
Our observations were compiled from four different sources: 1) Digital broad¬ 
band seismograms from the Saudi Arabian 1995/1996 PASSCAL temporary 
seismic network which include about 494 seismograms from earthquakes 
that occurred during the period December 31, 1995 to September 15, 1996 
and which represent about 50% of the Rayleigh-wave and about 60% of the 
Love-wave data; 2) Digital seismograms recorded by the permanent broad¬ 
band stations in the region for the period between 1990 and 1996 which rep¬ 
resent 26% of the Rayleigh-wave and 36% of the Love-wave data; 3) Analog 
observations of Rayleigh waves from the regional WWSSN stations recorded 
between 1970 and 1979 which represent 19% of the total Rayleigh-wave data; 
and 4) Analog observations from RYD station recorded between 1981 and 
1987 which represent 5% of both Rayleigh- and Love-wave data. 

In Figures 2 and 3 we present the great circle path coverage along which 
dispersion measurements were made for both Love- and Rayleigh-waves. A 
total of 987 Rayleigh-wave and 682 Love-wave paths were available. A maxi¬ 
mum of 916 rays for Rayleigh and 682 rays for Love were used in the period 
range of 11-13 seconds (this is the average of the total number of ray paths 
over the 3 periods). The number of great circle paths decreases for longer 
periods and reaches 376 for Rayleigh and 178 for Love at 56-60 s periods. 
Seismic stations used in this study are also shown in these Figures and 
include the XI95-96 IRIS PASSCAL digital broadband stations AFIF, HALM, 
RAYN, RIYD, RANI, TAIF, UQSK, and SODA; regional broadband GSN seis¬ 
mic stations ANTO, GNI, BGIO, KEG, and ATD; and the WWSSN stations 
from which analog data were digitized (TAB, SHI, EIL, JER, and RYD). 
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Figure 3. Distribution of ray path coverage for Rayleigh waves. 

The regional-scale tomographic maps presented below provide sig nifi cant 
constraints on the shear velocity and crustal thickness of the Arabian plate. 
These data were recorded mainly by seismic stations inside the plate or very 
close to its tectonic boundaries, and avoiding long paths that may traverse a 
number of different tectonic provinces and geological features. Surface waves 
propagating at regional distances are very useful in the period r ang e of inter¬ 
est in this study by contributing shorter periods and increased resolution. 
The stations deployed in Saudi Arabia provided significant coverage of the 
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Arabian shield. The data set also has significant coverage of the western part 
of Iran and all of Iraq. 

1.6 Methods of Analysis 

We present the results of regional tomographic inversion of dispersion 
data from Rayleigh and Love waves across Arabia using single station mea¬ 
surements of group velocity and applying the multiple filtering analysis tech¬ 
nique (Dziewonski et al., 1972, Herrmann, 1987). 

The dispersion measurements are obtained at period increments of 1 sec¬ 
ond between 5 and 20 seconds and at even periods only in the range 22-60 
seconds. Observations from each period are inverted separately and the 
images from adjacent periods were averaged. We parameterize the regional 
slowness variations using a uniform, one-degree by one-degree, grid of con¬ 
stant-slowness cells. Group velocity maps are produced using a conjugate- 
gradient least squares algorithm (Paige and Sanders, 1982). Laplacian 
smoothness constraints are incorporated in the inversion and thus we mini¬ 
mize a combination of group travel time misfit and a Laplacian measure of a 
two dimensional model roughness. The balance between group delay misfit 
and minimal roughness is selected empirically by running inversions with a 
range of smoothness importance weights and selecting the value that pro¬ 
duces the simplest model and still satisfactorily matches the observed group 
delays. The resulting velocity variations are presented as percent velocity 
perturbations from the average velocity of all the measurements (which was 
the initial model). 

We assume that the measured waves have propagated along the great 
circle path connecting the source and the receiver. Ritzwoller and Levshin 
(1998) discussed the problems expected to result from several factors such as 
off-great circle propagation, azimuthal anisotropy, and systematic event mis- 
locations near subducting slabs, and argued that these effects should not 
alter the tomographic maps of group velocities strongly beyond the resolution 
estimates. In addition, our use of short- distance paths helps minimize the 
likely deflection of the path from the great-circle arc. 

1.7 Resolution 

Estimating the resolution in a tomographic inversion is not a trivial task 
because resolution depends on complex factors such as the number of crossing 
rays, the density of sources and receivers, as well as random and possible sys¬ 
tematic uncertainty in the measurements. To estimate our resolving capabil¬ 
ity, we used standard "checker-board" tests. Specifically, we tested two mod¬ 
els. Each of the two models consists of 1° x 1° cells in which a velocity pertur¬ 
bation of +- 10% of the average was chosen to have a dimension of 3° x 3° 
for the first model and 8° x 8° for the second model (Figure 4). The checker¬ 
board test results show that features of dimension 8° x 8° are reasonably 
recovered, especially in the short period ranges where the number of ray 
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paths is maximum. In contrast, features of dimensions smaller than that are 
hardly recognized using the current coverage of the data available. The pre¬ 
cise amplitude of the anomaly is difficult to estimate due to dependence on 
damping and smoothing (which complicates tests involving models with 
sharp contrasts) but the pattern of variations is reasonably well recon¬ 
structed. We performed tests for both Love and Rayleigh waves and tested 
the resolution for the maximum and mi nimum number of ray paths in each 
case. 
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a) 


b) 


c) 


Figure 4. Checker board test results: a) The original velocity perturbation model using a fea¬ 
ture of 8° x 8° in dimension, b) Tomographic inversion results for Rayleigh-waves using the 
maximum number of rays, c) Rayleigh-waves results using the minimum number of rays. 


1.8 Group Velocity Variations in the Middle East 

In Figure 5 and Figure 6 we present the tomographic images constructed 
for Love and Rayleigh waves respectively. Each map shows the average 
group velocity variations for three adjacent periods and above each map we 
list the average group velocity for that period range and the average number 
of rays for the same periods. One striking observation in these maps is the 
consistency of the results from both Rayleigh and Love waves. The distribu¬ 
tion of faster and slower than average regions is strikingly similar in both 
data sets, especially in the short periods. It is clearly evident that the Ara¬ 
bian shield is characterized by relatively higher than average seismic veloc¬ 
ity, while the rest the Arabian platform, which is covered by an eastward- 
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thickening sedimentary section, has in general slower than average velocity. 
The boundary between the fast and slow regions is sharp and correlates well 
with the boundaries of the Arabian shield, especially at periods shorter than 
10 seconds. For longer periods, a fast region correlates well with the bound¬ 
aries of the Red Sea, especially for the Rayleigh waves. In general, the seis¬ 
mic velocity is higher in western Arabia than in the eastern or the northeast¬ 
ern parts of the plate. The resolution diminishes for long periods where the 
number of rays becomes smaller. The average group velocity of Love waves 
increases from about 2.44 km/s at 5-7 seconds to 3.98 km/s at 56-60 seconds, 
while that of Rayleigh waves increases from 2.38 km/s at 5-7 seconds to 3.74 
km/ s at 56-60 seconds. The mean velocities at our longer periods (50-60 s) are 
about 8% slow for Love waves, but only about 3% slow for Rayleigh waves. 

We do not have adequate coverage to image the southeastern part of the 
Arabian plate where it is covered mainly by the Empty Quarter region and 
the localized fast anomaly located east of the Gulf of Aden. The slow feature 
located in the eastern Mediterranean region should be interpreted as a result 
of bias in the data since there are not enough traverses crossing this region. 
The resolution of the oceanic structure south of the Arabian plate is limited 
to a small region that produces an apparent localized heterogeneity. Our 
results are consistent with those of Ghalib (1992) who analyzed analog seis¬ 
mograms recorded at stations TAB, SHI, EIL, and JER to outline the lateral 
variation of shear-wave velocity beneath the Arabian plate for depths from 5 
to 80 km. Ghalib (1992) concluded that shear wave velocity within the crust 
is higher in the shield region than in the platform area for depths less than 
40 km, but the pattern is reversed at depths below 40 km. This is consistent 
with the conclusion of Woodhouse and Dziewonski (1984) regarding the possi¬ 
ble existence of low-velocity anomalies in the upper mantle along the western 
and southwestern Arabian plate. 
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Figure 5. Tomography map of the Arabian plate Love-wave group velocity. The legend shows 
the % of the velocity perturbations. The average group velocity and the average number of 
rays used are given for each period interval. 
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Figure 6. Tomography map of the Arabian plate Rayleigh-wave group velocity. 
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2. CRITERIA FOR USING GENETIC ALGORITHM TO MODEL SUR¬ 
FACE-WAVE WAVEFORMS 

Tao-Ming Chang, Robert B. Herrmann, and Charles J. Amm on 

2.1 Abstract 

In this study, we evaluate the suitability of four criteria in modeling com¬ 
plicated waveforms such as surface waves. We design these four criteria 
using some direct measurable quantities such as Nl, N2 norm, normalized 
cross-correlation value at zero time lag, and time shift information of maxi¬ 
mum envelope and maximum amplitude of cross-correlation. The tests are 
performed using a genetic algorithm search method. Seventeen repeated 
experiments are conducted. Each experiment has four GA r uns with differ¬ 
ent criteria. The test results suggest that the criterion that combines phase 
and group velocity information is the best for modeling surface-wave wave¬ 
forms. The result of one randomly chosen experiment is shown and an over¬ 
simplified case is used to explain why this criterion is better than the others. 

2.2 Introduction 

Among seismological applications, the most common criterion of good- 
ness-of-fit is the L2 norm. This is because most seismological optimization 
problems are implemented using least squares minimization. For discrete 
data points, the least squares method and L2 norm are the natural choice. 
Indeed, for most studies the L2 norm behaves well and the resul ting surface- 
wave dispersion data inversion, modeling of teleseismic P waveforms, or 
receiver function inversion is acceptable. There is no need to further investi¬ 
gate the other possible criteria. 

Some people feel it is not totally convincing to use this simple criterion 
for every problem, especially when they try to use an L2 norm to model some¬ 
what complicated waveforms; for these cases they are starting to use more 
complicated criteria. For example, Zhu and Helmberger (1996) modified Zhao 
and Helmberger’s method (1994) to extract source information from broad¬ 
band regional seismograms by using combined LI and L2 norms. Other 
norms have been used to model the surface-wave waveforms. Lemer-Lam 
and Jordan (1983) and Cara and Leveque (1987) used cross-correlation tech¬ 
niques to formulate inversion schemes for inverting surface-wave waveforms. 
Therefore, we feel that the success of the simple L2 norm may be due to the 
simplicity of the problems, such as the use of travel times in determining 
hypocenter locations, teleseismic P waveforms in estimating source parame¬ 
ters, or receiver functions in inverting velocity structures underlying receiver 
sites. 

However, we believe there is a need to use complete waveforms to extract 
more information in order to better constrain earthquake source processes 
and propagation effects. Hence it is time to examine the suitability of 
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different criteria for the inverse problems. 

2.3 Purpose 

In this study, we will use the genetic algorithm (hereafter GA) as a tool 
to check the suitability of several criteria. It is not practical to test many cri¬ 
teria on different types of seismic signals, so we limit ourselves to modeling 
surface-wave waveforms. The test results will be useful only for this applica¬ 
tion. 

2.4 Basic Components For Designing Criteria of Goodness-of-Fit 

In implementing the GA search method, we note that the goodness-of-fit 
criterion is user defined, and can be a combination of a variety of different 
criteria. The design of a goodness-of-fit criterion for modeling a surface wave 
is very subjective and depends on the features of the surface-wave signal 
itself. To begin with the design of the criterion, it is necessary to understand 
the characteristics of the object which will be modeled. 

There are two major characteristics of surface waves. First, a surface- 
wave signal has a longer duration and a more complicated waveform than 
any single, pulse-like body-wave phase. When modeling such long-duration 
complicated waveforms, there is a cycle-skipping problem which may produce 
an unreasonably low or high velocity model. In addition, when processing 
surface-wave data, we cannot shift synthetic seismograms to match the 
observed arrival time, a well adopted technique in processing body wave data 
in the case of a receiver function inversion. Therefore, the cycle-skipping 
problem will be the first problem that will need to be resolved. Second, sur¬ 
face waves usually have a broad frequency content. Any single measurement 
from a waveform may not be adequate to represent the whole waveform 
because of the variability in amplitude levels as a function of frequency. For 
example, a single cross-correlation measurement only represents the similar¬ 
ity of the shape of the largest amplitudes, which are typically high-frequency 
for shallow crustal earthquakes in broadband data. This will only resolve the 
very shallow part of the structure and will leave the deeper structure uncer¬ 
tain with high variation. To overcome this problem, we can divide the fre¬ 
quency range of interest into several subranges and evaluate the goodness-of- 
fit of narrow-band filtered observed and synthetic seismograms for each sub¬ 
range. For our tests on regional seismograms, we divided the period range 
(10-50 sec) into 4 intervals : (10-20 sec), (20-30 sec), (30-40 sec), (40-50 sec). 
Using these period intervals as the ranges to bandpass filter observed and 
synthetic seismograms, a cross-correlation value is computed for each inter¬ 
val and a composite cross-correlation can be used as our goodness-of-fit. 

In this study, several different criteria have been tested as indicators of 
goodness-of-fit. Some of them are direct measurements from waveforms, the 
rest of them are a combination of those direct measurements. The direct 
measurements include the N1 norm (ATI), N2 norm (N 2), normalized cross- 
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correlation value at zero time lag (CC 0 ), the time shift of envelope of cross¬ 
correlation ( TS U ), and the time shift of maximum cross-correlation (TS C ) of 
windowed narrow-band filtered observed and synthetic seismograms. These 
basic criteria are briefly described. 

2.4.1 N1 Norm (Nl) 

For each sub-period range, the Nl norm is constructed from a normalized 
LI norm as 


_ Z W(tj) • 1 Fj* obsitj) - Fj * synitj) [ 

J EW(ti)- IFj^obsitJl {1) 

where W(£ £ ) is a cosine taper window which approximately has comers corre¬ 
sponding to group velocities of 4.8 and 2.6 km/sec, Fj is a band-pass filter, * 
means convolution, and obs(ti ) and syn(^) represent observed and synthetic 
seismograms respectively at time £;. The value from above equation will 
range from zero to some positive number. However in our GA search, the 
goodness-of-fit value is mapped into the range [0.0,1.0]. It is also necessary 
to map the Vj value into 1.0-0.0 range. The mapping procedure used in this 
study is as follows 


Nlj = 


l-log 10 (l + V>1.5 
HT 6 


0<Vj< 3.64158 
Vj > 3.64158 


( 2 ) 


Of course, the measurement can be performed for several wide-band frequen¬ 
cies (using a Butterworth filter) or several narrow-band frequencies (using a 
Gaussian filter). The goodness-of-fit value Nl is actually the combination of 
different frequency bands which emphasizes a solution that fits all frequency 
bands well. In the following, we use a geometric instead of arithmetic mean 
for goodness values of different frequency bands. 
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In a more complete form, the Nl norm used is 
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2.4.2 N2 Norm ( N2 ) 

Using the same mapping procedure as for the Nl norm, the N2 based 
norm is defined as 
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N 2 = 
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Like our N1 norm, this norm equally weights different sub-frequency bands 
instead of the domin an t frequency band, but internally uses an L2 norm. 


2.4.3 The Cross-Correlation Value at Zero Time Lag (CC 0 ) 

To model a surface-wave waveform, we hope that for every frequency 
range the normalized cross-correlation value at zero time lag equals unity, 
which means the synthetic seismogram is perfectly aligned with the observed 
seismogram. The normalized cross-correlation value ranges from -1 to 1. The 
following mapping procedure will emphasize the in-phase part. 


CC 0 = 


fl + <D(0)f 

2 J 


( 6 ) 


where O = [WitJFj * obsit^iWit^Fj * syn(ti )] is the normalized cross-corre¬ 
lation of windowed narrow-band filtered observed and synthetic seismo¬ 
grams. 

2.4.4 The Time Shift of Maximum Envelope Value of Cross-Correla¬ 
tion (TSJ 

The physical meaning of the time shift of the maximum envelope value of 
the cross-correlation of windowed narrow-band filtered observed and syn¬ 
thetic seismograms can be closely related to the difference of group velocity 
between observed and synthetic seismograms at that particular frequency. 
For modeling a surface-wave waveform, the best solution is that for all fre¬ 
quency subintervals there are no time shifts between observed and synthetic 
seismograms. The time shift (TS U .) for each sub-period interval is measured 
from the m aximum envelope value of the cross-correlation away from the zero 
time lag. The averaged absolute time shift will range from zero to some posi¬ 
tive number. Here another mapping procedure, different from the N2 norm 
mapping procedure, is needed to map time shift information into the good¬ 
ness range (0,1). The mapping procedure will be as follows: Assume surface 
waves travel with a group velocity of about 3.5 km/sec. First compute the 
percentage of time-shift with respect to the total travel time. Then, map it 
into the goodness range (0,1) using the following equation: 


TS U = 



f >1 
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(7) 


This criterion guarantees that for each subinterval the goodness value will be 
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equal to one when there is no time shift and the goodness will be 0.5 when 
the travel time difference reaches roughly 7%. The function TS U penalizes a 
poor fit. 

2.4.5 The Time Shift of Maximum Value of Cross-Correlation (TS C ) 

This criterion is closely related to phase velocity. Similar to TS U , when 
modeling a surface wave, we wish this time shift to be zero. The time shift is 
measured from the maximum value of the cross-correlation with respect to 
zero time lag. Using the same mapping procedure as TS U , the TS C is defined 
by the following equation. 



2.5 Four Combined Criteria 

Having defined several functionals that characterize a fit, we now create 
some combinations for testing overall fit. 

2.5.1 Criteria 1 

The first criterion which will be tested is (9). The concept behind the 
design of this criterion is that we want the GA search to focus on the models 
that minimize both LI and L2 based norms for every sub-frequency band. 
Initially, the N2 norm will be strongly affected by the biggest amplitude sig¬ 
nal, but the N1 norm will be less influenced by the largest amplitudes and 
will permit the small amplitude seismic signals to play a slightly more impor¬ 
tant role. 


l 

{N1 ■ N2}2 


(9) 


2.5.2 Criteria 2 

The second criterion is given by (10). It is well-known that using a sin¬ 
gle type of surface-wave dispersion data (such as group velocity or phase 
velocity) in an inversion will not necessarily match the waveform. This is 
because the inverted model which matches group velocity data may have a 
phase change, and results that match phase velocity data may still have a 
group delay. This is why we combined both types of information in this crite¬ 
rion. We hope the GA search will find models that have no time shifts in 
group and phase velocity for every sub-frequency interval. 

l 

f 1 2 


TS U ■ TS C 


( 10 ) 



2.5.3 Criteria 3 


The third criterion is (11). This criterion is similar to the second crite¬ 
rion but has a CC 0 term. The CC 0 contains redundant information because 
when a model has no time shift between observed and synthetic seismo¬ 
grams, the normalized cross-correlation value at zero time lag will be unity. 
The design of this criterion is for a comparison with Criteria 2 to see whether 
extra information will help the search converge faster. 
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f 1 2 
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2.5.4 Criteria 4 

The last criterion is (12). This is a combination of Criteria 1 and 3. We 
wish to consider the performance when a criterion utilizes all the possible 
information. 


| [N1 ■ N2]z 


CC 0 -(TS u TS c )l 



( 12 ) 


The criteria constructed have as their goals a characterization of a good 
visual effect and a sensitivity to sharply reject models that do not provide a 
good fit. 

2.6 Test Results 

It is impossible to reach a definitive conclusion based on one experiment. 
We use these four criteria in a GA search. We repeated this testing procedure 
for 17 experiments with small variations in the starting model. The number 
of experiments may not form a significant base from a statistical viewpoint, 
but we reach some understanding about the capabilities of different criteria. 
In the following section, we randomly choose one experiment result as an 
example to show the search results. 

The experiment has four GA search tests which use the different criteria. 
For each search test, we will show groups of three figures. The first figure in 
a test will show the ten best models, their waveform comparisons with the 
observed seismogram (filtered between 10 and 100 seconds and displayed 
from 50 to 250 seconds), and the distribution of searched models with good¬ 
ness values greater than a certain threshold. The second figure will show the 
seismograms of the ten best final models filtered in 4 period-ranges, (10-20 
sec), (20-30 sec), (30-40 sec), and (40-60 sec), displayed in a 100 to 300 second 
travel-time window together with the observed seismogram. The third figure 
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consists of two parts. The first part shows the group velocity extracted from 
observed and synthetic seismograms of the ten best final models u sing the 
multiple filter technique (Dziewonski et ah, 1969; Herrmann, 1973). The sec¬ 
ond part shows the averaged velocities of the top 30, 40, 60, and 80 km of all 
searched models with respect to goodness-of-fit values. 

The first GA search test used the Criteria 1. The results are shown in 
Figures 7, 8, and 9. In Figure 7, we can see the peak amplitude of the seis¬ 
mogram of the best model is time aligned better than the sixth best model. 
However from Figure 8, we can see that the first model does better in the 
20-30 second range and that this is the reason why the first model has a 
higher goodness value than the sixth model. Figure 9 shows two empirical 
ways to present the uncertainty of searched models. The group velocity curve 
can indicate the period range that deviates most from the observed curve. 
The shear velocities averaged over the upper 30, 40, 60, and 80 km show 
whether the GA searched models converge to a value which can be used to 
reconstruct the velocity profile. 

The second GA search test used the Criteria 2. The results are shown in 
Figures 10, 11, and 12. The synthetic seismogram of the best searched model 
fits the observed waveform very well with only a slightly slow higher mode. 
In Figure 12, we see that these ten models have their group velocity curves 
within ±0.1 km/sec range of the observed curve. However, for periods less 
than 10 seconds, none of the ten best models match the very low group veloc¬ 
ity. That can be explained because the "observed" model has a very low veloc¬ 
ity sediment layer at the top which may produce this very low group velocity 
for periods less than 10 seconds, which is beyond our modeling range. 

The third GA search test used the Criteria 3. The results are shown in 
Figures 13,14, and 15. 

The last GA search test used Criteria 4, with results shown in Figures 16 

-18. 

Although all these criteria can find models that fit the waveform reason¬ 
ably well, we need a way to compare the effectiveness of the models resulting 
from these four different criteria. Note that goodness-of-fit values the from 
different criteria are not directly comparable. To properly compare, we use 
the best model from each GA search test and evaluate it against each of the 
four criteria. The results are listed in Table 1. To understand this table, we 
would expect a criteria to give the best possible model, so that the diagonal 
would be dominant. In Table 1, the four evolutions of the best model of each 
GA search test are placed in one column. For a given model, the values of the 
other Criteria fill in the row. Table 1 indicates that the Criteria 3 model (row 
3) has the best value under TEST 2, 3 and 4 and is marginally poorer under 
TEST 1 than the model for Criteria 2 (2nd row). 
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Figure 7. The GA search results using Criteria 1. The waveforms of ten final best searched 
models (gray) are compared with the observed waveform (black). The displayed seismograms 
are filtered between 10 and 100 seconds, and the displayed range is 50 to 250 seconds. The 
right bottom plot shows these ten models and the best one is plotted as a heavy black line. 
The left bottom plot shows searched models with goodness values greater than 0.5. Lighter 
shades of gray imply a better fit. The outer bounds of the search are also indicated. 
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Figure 8. The observed and synthetic seismograms of ten models have been filtered in four 
different period ranges. The top left part is filtered for 10 to 20 seconds. The top right part is 
filtered for 20 to 30 seconds. The bottom left part is filtered for 30 to 40 seconds. The bottom 
right part is filtered for 40 to 60 seconds. The seismograms are displayed for 100 to 300 sec¬ 
onds. 


Table 1. The goodness values for 4 GA search tests 


Model 

from 

TEST 1 
(Criteria 1) 

TEST 2 
(Criteria 2) 

TEST 3 
(Criteria 3) 

TEST 4 
(Criteria 4) 

Criteria 1 
Criteria 2 
Criteria 3 
Criteria 4 

0.851907 

0.974151 

0.972166 

0.910052 

0.873815 

0.961544 

0.971157 

0.921201 

0.899356 

0.978459 

0.981818 

0.939683 

0.860601 

0.945162 

0.953598 

0.905907 


To compare the abilities of criteria we introduce two voting schemes. In the 
first, we examine each row to seek the biggest value for that criteria and give 
a vote to that GA Search test. In the second, we check the diagonal values to 
see whether its value is the largest in that row. If answer is yes, we give a 
vote to that GA search test at that column. For example, from Table 1, TEST 


22 






© 

o 



o 



Figure 9. The top panel shows the group velocity (km/sec) information extracted using the 
MFT technique compared to dispersion curves from the best models. The observed value is 
shown as a black line. The bottom panel shows the shear velocity averaged over the upper 
30 (x), 40 (+), 60 (triangle), and 80 km (circle) of all GA searched models. 


3 which uses Criteria 3 in its GA search gets 4 votes for scheme 1 , and 1 
point for scheme 2; the other 3 tests get 0 votes under either scheme. We now 
restart the GA search process for 17 initial models for each of the four 
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Figure 10. GA search results using Criteria 2. See Figure 7 for explanation. 


goodness of fit criteria. The composite tally of votes is given in Table 2. 
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In Table 2, we can see that Criteria 2 and 3 are roughly equivalent and 
are better than Criteria 1 and 4 in determining a good model. 

Why does Criteria 2 seems to be better than Criteria 1? It must be due 
to reasons that force the hyper-surface of the object function to be smoother 
than the others. Therefore, the GA evolution mechanism will gradually 
approach the best goodness region without much trouble. Although it is 
impossible to show the object function surface because it is a multi-dimen¬ 
sional surface, it is also impossible to evaluate such a surface for the infinite 
number of models that need to be computed. All we can say is that the object 
function is smoother in some over-simplified cases that we can compute. For 
example, if observed and synthetic seismograms are very similar to each 
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Figure 12. Uncertainty for the second GA search test. See Figure 9 for explanation. 

other with only some time-shift between them, we can see the object function 
of Criteria 2 is smoother than the object function of Criteria 1. 

To illustrate this, we use a synthetic seismogram in Figure 19 for our 
’observed" data and shift this seismogram for our "synthetic" data and evalu¬ 
ate the goodness values for four different criteria. The bottom panel shows 
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narrow range of high goodness values. This is mostly controlled by the peri¬ 
odic behavior of the LI, L2 based norms and the cross-correlation values. 
However, Criteria 2, which only uses the group and phase velocity informa¬ 
tion, has a broader and smoother range of goodness values which means it is 
more tolerant of time shifts than the other three criteria. 

2.7 Conclusion 

In this study, we designed four criteria for modeling surface-wave wave¬ 
forms. These four criteria were designed using some directly measurable 
quantities such as LI and L2 norms, normalized cross-correlation value at 
zero time lag, and time shift information of maximum envelope and maxi¬ 
mum amplitude of cross-correlation. These four criteria were tested using 
the genetic algorithm search method for 17 repeated experiments. Each 
experiment consisted of four GA runs with different criteria. 

The results suggest that the second criteria, which only uses group and 
phase velocity information, is the best for modeling surface-wave waveforms. 
The other three criteria may be influenced by the periodic behavior of the LI 
and L2 norms and cross-correlation value and have smaller and narrower 
regions of high goodness values. This means the second criteria is more 


28 




o 

o 



o 



Figure 15. Uncertainty for the third GA search test. See Figure 9 for explanation. 

tolerant of incorrect models and can still use information of those models in 
the evolution process of genetic algorithm. 

Although we only performed 17 experiments, and this may not be enough 
from a statistical viewpoint, we believe the result of this study is an initial 
step m searching for a good criteria for use in modeling complicated seismic 
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Figure 17. Seismogram comparison for ten best GA searched models using criteria 4. See 
Figure 8 for explanation. 
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Figure 18. Uncertainty for fourth GA search test. See Figure 9 for explanation. 
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Figure 19. Illustration of object functions for four different criteria in an over-simplified case. 
See text for detail. 
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3. MODELING SURFACE-WAVE WAVEFORMS USING A GENETIC 

ALGORITHM 

Tao-Ming Chang and Robert B. Herrmann 

3.1 Abstract 

We apply a genetic algorithm search method to model regional surface- 
wave waveforms. In this application, we use some seismological constraints 
to regulate the test models which improve the use for small population sizes 
(20 individuals) and short numbers of generations (50 generations). The 
waveform goodness-of-fit is based on measurements of travel time differences 
of normalized cross-correlation of observed and synthetic seismograms that 
have been narrow-band filtered in subranges of the period interval of inter¬ 
est. The synthetic tests show our application of the genetic algorithm search 
method can successfully find several good final models that can be used as 
input for further detailed studies. The real data tests have been applied to 
data from the 1995 West Texas earthquake and have shown a satisfactory 
waveform fit. We have not found a good way to quantify the uncertainty of 
the genetic algorithm search method, but provide an indirect way to repre¬ 
sent the uncertainty qualitatively. 

3.2 Introduction 

Most seismological inverse problems are nonlinear. The techniques used 
to solve such nonlinear problems can be placed into two groups. The strategy 
of the first is to linearize nonlinear problems, then use iterative processes to 
seek a better solution by using the gradient information of the hypersurface 
of the misfit function. The strategy of the second is to directly search the 
model space and find the acceptable models. 

Methods such as least squares, steepest descents, and conjugate gradient 
belong to the first group. Although these methods are widely used in seismol¬ 
ogy, the requirement of a good starting model is a well-known disadvantage. 
For studying large scale, long wavelength features or deep features beneath 
the lithosphere, this would not be a real problem because the research of the 
past half century already provides some good starting models, e.g., PREM 
(Dziewonski and Anderson, 1981) and IASPEI91 (Kennett and Engdahl, 
1991). On the other hand, for studying crustal or lithospheric structure, a 
good starting model may not be available because the crust or lithosphere is 
the most structurally heterogeneous region in the Earth. We must investi¬ 
gate the structure of the crust and lithosphere all over the world to study its 
evolutionary history, to understand the tectonic processes, to better predict 
seismic activity, and to lower the damage of earthquake hazards. 

Therefore, it is necessary to find a way to investigate many possible mod¬ 
els, that is, to search the whole model space and to select the acceptable mod¬ 
els of structure and their variations. The Monte Carlo, simulated annealing 
(SA), and genetic algorithms (GA) search methods belong to this group. The 
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Monte Carlo method is a random search method that has been used in seis¬ 
mology for a long time (for example, by Keilis-Borok and Yanovskaja, 1967; 
Press, 1968). Press (1968) showed a successful experiment that used the 
Monte Carlo method to search for the model that can produce correct body- 
wave travel-times, surface-wave dispersion, the earth’s free oscillation peri¬ 
ods, mass, and moment of inertia. Six models were found from about five mil¬ 
lion randomly generated models. Examining these six models, we find that 
the structures for the lower mantle are consistent, but that large fluctuations 
occur in the upper mantle (lithosphere). As pointed out by Press (1968), the 
reason for using the Monte Carlo method is that it offers the advantage of 
exploring the range of possible solutions and indicates the degree of unique¬ 
ness achievable with currently available geophysical data. However, the 
Monte Carlo method has its own disadvantages. As discussed by Keilis- 
Borok and Yanovskaja (1967), the Monte Carlo method does not use informa¬ 
tion obtained from previous trials in the next trial. 

Recently, two better search algorithms which utilize the information of 
previous trials, simulated annealing and genetic algorithms, have become 
very popular in seismology. The simulated annealing method mimics the 
crystallizing process observed in chemistry. The genetic algorithm is inspired 
by the evolution process observed in biological science. The seismological 
applications have been the estimation of residual statics (Rothman, 1985, 
1986; Wilson and Vasudevan, 1991), waveform inversion of reflection data 
(Sen and Stoffa, 1991, 1992; Stoffa and Sen, 1991; Sambridge and Drijkonin- 
gen, 1992; Drijkoningen and White, 1995), earthquake hypocenter location 
determination (Sambridge and Gallagher, 1993), receiver function inversion 
(Shibutani, Sambridge, and Kennett, 1996; Zhao et al., 1996; Nerves et al., 
1996), determination of earthquake source parameters (Hartzell and Liu, 
1995; Zeng and Anderson, 1996; Kobayashi and Nakanishi, 1994), estimation 
of the source time function (Courboulex et al., 1996). Simulated annealing 
and genetic algorithms have also been used to study many different type of 
data such as teleseismic P waves (Steck, 1995; Zhou et al., 1995), surface- 
wave dispersion data (Lomax and Snieder, 1994, 1995; Yamanaka and Ishida, 
1996), and seismic exploration data (Jervis et al., 1993, 1996; Jin and 
Madariaga, 1994; Mallick, 1995; Nolte and Frazer, 1994), and other applica¬ 
tions (Pullammanappallil and Louie, 1994; Riepl et al., 1995; Vasudevan et 
al.; Velis and Ulrych, 1996). There is one thing in common in all these appli¬ 
cations: they all deal with nonlinear problems (Scales et al., 1992; Gallagher 
et al., 1991). In addition to their wide range of applications in seismology, GA 
and SA may also provide an intuitive way to analyze the error or uncertainty. 
For example, in receiver function inversion, Ammon et al. (1990) showed that 
the final models were dependent on the initial models. Shibutani et al. (1996) 
showed that genetic algorithm can estimate an average model which is more 
stable and less dependent on the starting model. 
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We see the goal of the global search techniques is to obtain several good 
starting models and the uncertainty range, which will then be refined using 
using iterative waveform inversion techniques. In this paper, we will focus on 
several topics to explain why we study the surface wave, why we use the 
genetic algori thm instead of simulated annealing, and what our purposes 
and expectations are. 

3.3 Purpose 

At present the Genetic Algorithm has been applied in many topics in 
seismology such as hypocenter location determination, receiver function 
inversion, est im ation of the source time function and surface wave dispersion 
data inversion. But, the results show there is no guarantee of reaching global 
minimum by using GA. So, the proper questions about the use of GA in seis¬ 
mology are: What is the strength of genetic algorithm method? What can be 
accomplished by using GA? How do we utilize this tool? 

In this study, we use GA as a tool to quickly delineate the possible mod¬ 
els for surface-wave waveform modeling. Under some reasonable a priori 
conditions, we use GA to search the model space and to determine a few good 
starting models for other waveform inversion algorithms which are better for 
such a task. 

Because the genetic algorithm crudely searches the model space, the 
search results can be used to show the possible ranges of uncertainty. This is 
important because most goodness-of-fit measures used in seismology are not 
useful in indicating the actual uncertainties. For example, in damped least 
square inversion, the resolving kernel and standard deviation are dependent 
on the damping value. It is possible to adjust the damping value and show a 
different range of uncertainty. Therefore, it would be better to have the 
search results indicate the possible uncertainty. In this study, we plot the 
search results to show the deviation range of models. Although we can only 
provide qualitative instead of quantitative information, such presentation 
will provide intuitive information in the manner of the "checker board" test 
used in tomographic inversion. Also, as pointed out by Sambridge and Dri- 
jkoningen (1992), the development of error analysis for a GA search remains 
a great challenge if one is only to rely on GA. 

3.4 Why Study Surface Waves? 

As stated previously, the main goal is to establish a set of algorithms for 
studying lithospheric structure. Although different seismic phases can be 
used for such a task, we select surface waves for several major reasons. First 
of all, high quality surface-wave data are easily obtainable because of widely 
deployed broadband instruments. For most shallow earthquakes, the surface 
wave usually has a higher signal/noise ratio than the body-wave signal. Sec¬ 
ond, for any local small aperture broadband seismic network, the records may 
form a huge database for studying the regional structure. However, the 
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difficulty in such research is that signals may have a higher frequency con¬ 
tent than signals used in global tomography studies. At these higher fre¬ 
quencies, the surface wave travels mostly through the crust which is the 
most heterogeneous part of the earth. It is really necessary to have a method 
to obtain reasonable starting models for detailed waveform modeling. 
Finally, since surface waves sample the near surface structure, modeling sur¬ 
face waves using a 1-D model can provide an average structure. Although this 
does not sound exciting, it should serve as the foundation for studying compli¬ 
cated 3-D structures. As an analogy, seismic exploration, researchers are try¬ 
ing to extract background velocity (1-D), and use it as a starting point for 
studying more complicated structures (Koren et al., 1991; Bunks et al., 1995; 
Lauda et al., 1989). 

Here, we also would like to explain why are we using surface-wave wave¬ 
form data instead of dispersion data. Theoretically, if the data have a high 
S/N ratio, the waveform and dispersion data should be equally suitable for 
inversion. There are some reasons for us to use the waveform data. First, 
there has been a tremendous improvement in computation ability. Calculat¬ 
ing synthetic seismograms may only take seconds instead of the hours previ¬ 
ously. Second, the dispersion data were usually obtained using the multiple 
filter analysis technique. Sometimes unavoidable measuring errors will be 
introduced in such measuring procedures. In order to get the phase velocity 
data, another measuring process such as a phase matched filter is necessary. 
Of course, some artifacts may be introduced. It is very difficult to identify the 
dispersion for higher modes. We will try to find a good way to extract infor¬ 
mation directly from the waveform to avoid these problems. 

3.5 Comparison of SA and GA 

From the variety of applications in earthquake and exploration seismol¬ 
ogy, we find that both simulated annealing (SA) and genetic algorithms (GA) 
are good ways to perform the uncertainty assessment in a complicated non¬ 
linear problem. The algorithm chosen depends on the characteristics of prob¬ 
lem. As stated by Sambridge and Drijkoningen (1992) concerning the SA and 
GA methods: "any problem feasible by one could also be tackled by the other". 
To distinguish which of the two methods (GA or SA) is better for the surface- 
wave waveform modeling problem, we need to understand both algorithms 
and the purpose of our application. As is well known, generating multi-mode 
surface-wave synthetics is computationally intensive. Therefore, the compu¬ 
tation time will be a crucial factor for selection of the algorithm. Of course, 
there are some applications of an artificial neural network algorithm in seis¬ 
mology (Leach et al., 1993; McCormack et al., 1993; Paulli and Dysart, 1990; 
Wang and Teng, 1995, 1997; Wang and Mendel, 1992) but they are less popu¬ 
lar than GA and SA and may require very large training sets. 
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3.5.1 Workflow of Simulated Annealing Method 


The computational procedure of simulated annealing consists of the fol¬ 
lowing steps: 

• start from an arbitrary model 


• temperature-loop : at temperature T - T 0 -k ■ ST 

□ parameter-loop : for model parameter i = 1, • • •, m 

fix all other parameter value except 

for the parameter S t , there are n possible values. 

o possible-value-loop : Sy, j = 1, • • •, n 

calculate an energy function E(Sy) (the normalized cross¬ 
correlation of observed and synthetic seismogram is one 
example) 

o end of possible-value-loop 

□ end of parameter-loop 
compute the probability distribution 


W,) = - 


r E(Sy) 
ex Pt f ~ 1 


X ex Pf—sr^-1 

i=i 1 


(13) 


• end of temperature-loop 

We can see that for each temperature, it is necessary to perform (m • n) for¬ 
ward computations of synthetics. If there are k steps in lowering tempera¬ 
ture to reach the global minimum, the total forward computation will be 
(k • m • n). The problem is that there is no rule for choosing the starting tem¬ 
perature T 0 and increment of temperature difference ST. Basu and Frazer 
(1990) designed a sequence of test runs to find the critical temperature. In 
spite of this, it is still too time-consuming for surface-wave waveform model¬ 
ing. 

3.5.2 Workflow of the Genetic Algorithm 

For an m member society evolving through n generation, the computa¬ 
tional sequence of the genetic algorithm is as follows: 

• Randomly generate m individual models as the first generation. 

• Generation-loop : for generation = 1, • • •, n 

Compute m synthetic seismograms (individuals) 

Evaluate each individual’s performance; i.e. calculate the goodness- 
of-fit. 

□ population-loop : for child = 1, • • •, m 

Based on the individuals’ performance (probability) 
select them as parents; 
change parents’ DNA; and 
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apply possible mutation on their children. 

□ end of population-loop 

• end of generation-loop 

We see it is possible to perform GA on a small population. This will be com¬ 
putationally more efficient than SA in the surface-wave waveform modeling 
problem. However, performing GA on a small population society has its own 
risk, as pointed out by Sambridge and Drijkoningen (1992). If the society’s 
members are not close to the global solution, the few relatively good individu¬ 
als in the society will multiply themselves and dominate the population (that 
is, the solution will be trapped in a local minimum leading to premature 
selection). This problem can be addressed by increasing the population size. 
However, our purpose for using GA in surface-wave waveform modeling is not 
to rely on the GA to reach the global minimum. Instead, we prefer to have 
several runs to see the possible uncertainty, get a rough idea about the struc¬ 
ture, and find some good initial models for other inversion techniques. In our 
test, a GA run will take 1 to 3 hours of CPU time on a 167 MHz workstation. 
It is affordable to have several reruns if we find it trapped in a local mini¬ 
mum. 

We thus select the genetic algorithm for modeling surface-wave wave¬ 
forms. 

3.6 Technical 

The basic concepts of the genetic algorithm are well presented by other 
authors (e.g. Sambridge and Drijkoningen, 1992; Stoffa and Sen, 1991). Here 
we will only address the technical issues in implementing GA to surface-wave 
waveform modeling. 

3.6.1 Generation Number and Population Size 

Although GA is a global search method which can potentially find the 
global minimum , we did not set that as the goal in this study. Due to the 
intensive computational load of generating multi-mode surface-wave synthet¬ 
ics, we limited our computations to a small population size and propagated it 
through only a finite number of generations. We hoped, by using the GA 
search method, to get some starting models for other inversion algorithms. 
To un derstand what generation number is sufficient for our purpose, we have 
tested the consequences of a large generation number. In this test, with 
results shown in Figure 20a, we propagated 500 generations and find that 
after 50 generations model improvement is less rapid, indicating a degree of 
convergence. Therefore, in the subsequent tests, we will only use 50 genera¬ 
tions for a small population (20). This will only consume 1 to 3 hours of CPU 
time on a 167 MHz workstation. In Figure 20b, we perform four repeated 
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Figure 20. (a) Top panel. This experiment shows the best fitness found by a small popula¬ 
tion (20 individuals) GA search through 500 generations. The result shows that 50 genera¬ 
tions is acceptable for the reason of saving computation time, (b) bottom panel. Comparison 
of 4 separate GA searches. Again, the 50 generation limit seems acceptable. 


experiments and. show their best fitness values through 50 generations. The 
experiments show that the GA search method is reliable and repeated experi¬ 
ments show similar evolution trends. 
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3.6.2 Smoothing Mechanism or Model Generating Mechanism 

The GA searches the model space randomly, therefore it has an advan¬ 
tage that a broad range of models can be tested and it is possible to reach the 
global minimum if the problem has been properly formulated and the search 
process has lasted long enough. However, because of its random search 
behavior, many unreasonable models can be generated which are not accept¬ 
able under seismological constraints. In this study, we found it necessary to 
use some seismological constraints to regulate the randomly generated mod¬ 
els. We call this process "smoothing". This process will help the GA find 
many "reasonable" models in such a short run for the small population size 
described previously. This smoothing concept is not new to this report; 
rather, several papers already mention similar ideas. For example, Lomax 
and Snieder (1995) used a "centrally weighted five-point" smoothing mecha¬ 
nism to smooth over their 18-node velocity profiles. Another example is from 
Sambridge and Drijkoningen (1992), who allow only models which have posi¬ 
tive velocity gradients to be generated when they modeled marine seismic 
refraction waveforms. Although these smoothing mechanisms depend on 
the parametrization, we believe this concept will be very important for seis¬ 
mologists who wish to apply GA to other seismological problems. Therefore 
we will briefly describe the smoothing process. First, we need to define a 
"reasonable model". For a regional shallow event, the surface wave can prob¬ 
ably only resolve the crustal structure. When low-pass-filtered the surface 
wave signal is a simple pulse. We will not see very detailed structure by mod¬ 
eling such signals. Under such limitations, we can only hope to find some 
velocity models that vary simply with depth. For this study, this is our defini¬ 
tion of "reasonable model". 

As shown in Figure 21, the free parameters will randomly generate 
models with inherent "zig-zag" behavior. Because of the period range of sur¬ 
face wave and the scale of structure, we wish to use a smooth velocity-depth 
curve instead of randomly generated models with a "zig-zag" pattern. There¬ 
fore, we smooth it in two ways that will be chosen randomly. The first way is 
to reduce the velocity contrast between any two adjacent layers without 
changing the total vertical travel time. The second way is to use a center- 
weighted triangular shape moving over the velocity profile to smooth the 
velocity model. This is similar to the smoothing mechanism used by Lomax 
and Snieder (1995). 

A constraint was used to check the validity of such a model. From most 
published models, we conclude that it is valid that the average velocity of the 
upper crust is lower than the average velocity of lower crust, that the average 
upper mantle velocity is faster than the average lower crust velocity, and that 
the upper mantle does not have strong velocity gradients. For example, if the 
Moho depth is 40 km, we will divide the first 100 km into four layers: layer 1 
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Figure 21. This figure shows how the smoothing mechanisms smooth the randomly gener¬ 
ated model (solid line) which always has a strong "zig-zag” pattern. The long dashed line 
shows the smoothed result by reducing velocity contrast. The short dashed line shows the 
smoothed result by moving a center-weighted triangular weighting function over the model. 

(0-20 km; upper crust), layer 2 (20-40 km; lower crust), layer 3 (40-60 km; 
upper mantle), layer 4 (60-100 km). The average velocity for the first three 
layers should increase monotonically, and the average velocity of the fourth 
layer should not deviate too much from the third layer (typically, it should be 
within 0.15 km/sec). 

As mentioned previously, the surface-wave period range generated by a 
regional shallow event may provide very little resolution for structure deeper 
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than the Moho. However, it is necessary to provide the deeper structure to 
fit higher modes or for modeling teleseismic seismograms. For this part of 
the velocity model, we do not generate it randomly. Instead, we construct a 
sequence of models (Figure 22) based on the SNA, TNA (Grand and Helm- 
berger, 1984), and IASP91 models (Kennett and Engdahl, 1991). The deeper 
model selected will be based on the average velocity of the fourth layer (usu¬ 
ally 60-100 km) mentioned previously. 


3.6.3 The Criteria of Goodness-of-Fit 

In GA, the goodness-of-fit criteria is user defined. It is a subjective crite¬ 
ria of goodness-of-fit for modeling surface waves. To start the design of the 
criteria, it is necessary to understand the characteristics of the object that 
will be modeled. 

There are two major characteristics of a surface wave. First, a surface- 
wave signal has a longer duration and a more complicated waveform behavior 
than any single, pulse-like body-wave phase. When modeling such long-dura¬ 
tion complicated waveforms, there is a cycle-skipping problem which may pro¬ 
duce an unreasonably low or high velocity model. In addition, when process¬ 
ing surface-wave data, we cannot shift synthetic seismograms to match the 
observed arrival time, a well adapted technique in receiver function inversion 
of body wave data. Therefore, the cycle-skipping problem will be the first 
problem to be resolved. Second, surface-waves usually have a broad fre¬ 
quency content. Any single measurement from a waveform may not be ade¬ 
quate to represent the whole waveform because of the frequency dependence 
of the signal content. For example, a single cross-correlation measurement 
only represents the similarity of the shape of the largest amplitudes, which 
are typically high-frequency for shallow crustal earthquakes. This will only 
resolve the very shallow part of the structure but leave the deeper structure 
uncertain with high variation. 

To overcome these two problems, we design our criteria of waveform 
goodness-of-fit as 
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First we divided period range of interest into several subranges and narrow 
bandpass-filter observed and synthetic seismogr ams for each subrange. For 
example, for a regional seismogram, we divided the period range (10-50 sec) 
into 4 intervals : (10-20 sec), (20-30 sec), (30-40 sec), and (40-50 sec). For each 
subrange period interval, we compute the cross-correlation of these narrow- 
band filtered observed and synthetic seismograms. After that, two important 
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Figure 22. Because the surface-wave generated by shallow focus earthquakes usually do not 
generate much long-period signal that can be used to resolve the deep structure, it is neces¬ 
sary to ensure that the deep part of the model be able to model teleseismic seismograms. We 
construct a sequence of possible models which are based on SNA, TNA, and IASPEI91 mod¬ 
els for the GA search method. 


time shift parameters can be measured: the time shift of the envelope of 
cross-correlation ( TS U ) and the time shift of maximum cross-correlation 
(TS C ). These two measurements are similar to surface-wave group and phase 
velocity differences between observed and synthetic seismograms. Assuming 
that the averaged surface-wave traveling speed is 3.5 km/sec and using epi- 
central distance ( dist ), we can roughly estimate the percentage of travel time 
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difference between observed and synthetic seismograms. Finally, the value of 
goodness-of-fit can be calculated. When TS C and TS U approach zero for each 
subrange, we can assume that a good solution has been found and its good¬ 
ness-of-fit is close to one. Otherwise, the exponential term in (13) will pull 
the goodness-of-fit away from unity for just a few percent travel time differ¬ 
ence. This criteria will also tolerate an incorrect Q model used to generate a 
synthetic seismogram because it is based on the normalized cross-correlation 
measurements. 

3.7 Synthetic Test 

We conducted several synthetic tests to examine the suitability of the GA 
search method. Firstly, correct and incorrect crustal Q models have been 
tested. In our experiment, we set averaged crustal Q s = 800, which deviated 
100% from the correct Q s = 400. As discussed in the previous section, the cri¬ 
teria of goodness-of-fit is based on normalized cross-correlation measure¬ 
ments, therefore it won’t show too much difference for the two cases. Sec¬ 
ondly, we have tested the consequence of using Butterworth and Gaussian 
filters. The results show roughly the same results for both filters. From 
these two experiments, we conclude that the design of goodness-of-fit criteria 
is very robust. The success of the GA search method depends only on the 
mechanism of the genetic algorithm (that is, selection, exchange DNA, and 
mutation) or the mechanism of generating models (that is, smoothing, param¬ 
eterization). 

We also conducted a test to see whether the GA search method has the 
ability to distinguish between a gradient Moho and a sharp Moho. Figures 23 
and 24 show the GA searched results for synthetic gradient Moho and sharp 
Moho. These tests assume the receiver is 500 km away from the source and 
that we have imperfect knowledge about crustal Q models. Each figure has 
three small plots. The top panel shows the ten best search results (gray line) 
compared with "observed" seismogram (black line) with their fitness value on 
the right; the displayed time interval is from 50 to 250 seconds. The bottom 
right plot shows these ten searched models (light gray line), the best model 
(solid black line), the "true" model (black dashed line), and search bound 
(thick dark gray line). Although, these ten models have somewhat good 
waveform fits when judged visually, their models span a certain width on 
each layer which may be used as an indicator for uncertainty. For example, 
the uncertainty around the Moho is large in comparison to other depths. This 
indicates that the surface wave has very little ability to resolve Moho struc¬ 
ture for such a short distance (500 km) and period range (10-100 sec). Also, 
the top layer of very low velocity sediment is not resolvable from this period 
range. The bottom left plots of Figures 23 and 24 show those search models 
that have their goodness-of-fit value greater than a certain threshold. The 
color indicating different goodness-of-fit values, for example very dark gray 
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for those model greater equal 0.9, medium gray for those between 0.7 and 0.9, 
gray for those between 0.5 and 0.7, and light gray for those below 0.5. From 
these experiments, we conclude that the GA search method can roughly 
resolve the shear velocity structure, but it needs more a priori knowledge to 
resolve the top sediment layer and Moho structure. 
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Figure 23. This figure shows GA test results for a gradient Moho model. The top panel 
shows waveform comparison (gray) for ten best searched models with their goodness-of-fit 
values on the right, against the original waveform (black). The bottom right plot shows the 
search bounds (outer gray), the ,, true ,, model (black dashed), ten best models (light gray), and 
the best model (black). The bottom left plot shows most models with lighter color, indicating 
a better goodness-of-fit value. See text for details. 
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Figure 24. This figure shows GA test results for a sharp Moho model. See Figure 23 for 
explanations and text for details. 


3.8 Real Data Test (1995 Texas Earthquake) 

We wish now to apply the GA search method to real data. We selected 
the April 14, 1995 western Texas earthquake to test. Here we will show the 
search results for two stations (WMOK and FCC) which have very different 
epicentral distances, 655 and 3244 km, respectively. Figure 25 shows the GA 
search results for WMOK. The displayed seismograms are filtered between 
10-100 seconds and cut between 100-300 seconds after origin time. From the 
bottom right plot, we see that the search results have a large uncertainty for 
structure below 40 km and a slightly lower velocity for the deeper part. This 
is reasonable because the observed data do not have enough long-period 
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signal to resolve the deep structures. However, if we focus on the shallow 
part of the models, we see well-resolved shear velocity structure for the top 
30 km which is contributed by the strong fundamental mode. The flaw of this 
search result is that none of the best ten models can fit the higher modes 
well. The GA model solution may be used as a starting point for other inver¬ 
sion techniques. 
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Figure 25. This figure shows GA search results for WMOK. See Figure 23 for explanations 
and text for details. 


Figure 26 shows the GA search results for FCC. The displayed seismo¬ 
grams are filtered between 20-100 seconds and cut to 750-1300 seconds after 
origin time. In this search, we divided the period interval of interest (10-100 
sec) into four sub-intervals: (10-20 sec), (20-35 sec), (35-60 sec), and (60-100 
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sec). The search results do not fit the observed waveform perfectly. The best 
model has the correct velocity for the low frequency signals which means the 
deep structure of its model is close to the truth but it fails to produce the cor¬ 
rect amplitude around the Airy phase. The tenth model has the correct 
amplitude and envelope shape around the Airy phase but fails to fit the long- 
period part. Therefore, we see all ten best search results could be used as a 
starting point for more detailed studies of waveform behavior. 
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Figure 26. This figure show GA searched results for FCC. See Figure 23 for explanations 
and text for details. 
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3.9 Discussion 

This study has shown that the GA search method is capable of finding 
several models that can fit the surface-wave waveform reasonably well. The 
GA search method can be applied to regional and teleseismic records because 
we construct a sequence of models for the deep structures. This method is 
very fast and therefore is practical for application to real earthquake data. 
The search models can be used as starting models for other detailed studies. 
The research goal has been fulfilled successfully. Although there is no quanti¬ 
tative method for defining the uncertainty of searched results, several sugges¬ 
tions have been made for indirect measures. 

What is the weakness of this tool? It is obvious that not all of the 
searches will be successful. Although we believe our criteria of goodness-of-fit 
is very tolerant of an incorrect Q model and is reasonably well behaved, it 
definitely is not the perfect one. That is why the GA search will sometimes 
yield excellent results and sometimes not. In this study, we feel that the find¬ 
ing of final good solutions is more dependent on the model-generating mecha¬ 
nism than on the mechanism of the genetic algorithm (evolution process). In 
other words, the more seismological constraints that have been applied, the 
greater efficiency and accuracy of the process. Also, if we know more about 
waveform behavior, we may design a more intelligent and efficient guided 
search method. 
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